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ABSTRACT 

We present results from Spitzer Space Telescope observations of the mid-infrared phase 
variations of three short-period extrasolar planetary systems: HD 209458, HD 179949 
and 51 Peg. We gathered Infrared Array Camera (IRAC) images in multiple wavebands 
(3.6 /im or 4.5 /im, and 8 /Ltm) at eight phases of each planet's orbit. We find the 
uncertainty in relative photometry from one epoch to the next to be significantly 
larger than the photon counting error at 3.6 /^m and 4.5 ^m. We are able to place 2a 
upper limits of only ~ 2% on the phase variations at these wavelengths. At 8 /im the 
epoch-to-epoch systematic uncertainty is comparable to the photon counting noise 
and we detect a phase function for HD 179949 which is in phase with the planet's 
orbit and with a relative peak-to-trough amplitude of 0.00141(33). Assuming that 
HD 179949b has a radius Rj < R p < 1.2Rj, it must recirculate less than 21% of 
incident stellar energy to its night side at the la level (less than 26% at the 2a level, 
where 50% signifies full recirculation). If the planet has a small Bond albedo, it must 
have a mass less than 2.4Mj (la). We do not detect phase variations for the other 
two systems but we do place the following 2a upper limits: 0.0007 for 51 Peg, and 
0.0015 for HD 209458. Due to its edge-on configuration, the upper limit for HD 209458 
translates, with appropriate assumptions about Bond albedo, into a lower limit on the 
recirculation occuring in the planet's atmosphere. HD 209458b must recirculate at 
least 32% of incident stellar energy to its night side, at the la level (at least 16% at 
the 2a level) , which is consistent with other constraints on recirculation from the depth 
of secondary eclipse depth at 8 /xm and the low optical albedo. These data indicate 
that different Hot Jupiter planets may experience different recirculation efficiencies. 

Key words: methods: data analysis - stars: planetary systems - techniques: photo- 
metric 



BACKGROUND 



1.1 Introduction 



Since the discovery o f the first Hot Jupiter system 
l|Mavor &: Quelozl Il995h and the subsequent realiz ation 
that such systems are ubiquitous (jButler et alj 120061 ) , as- 
tronomers have struggled to observe these planets and to 
compare these observations to theoretical predictions about 
their surface temperature variations, composition, albedo, 
and internal structure. The observation of extrasolar plane- 
tary t ransits, starting wit h HD 209458b (|Charbonneau et al.l 
120001 : iHenrv et al.l [§000) , have allowed for the determina- 
tion of the mass and size of these planets, but it was 
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only with the observat ion of secondary e clipses, start- 
ing with HD 209458b llDeming et alj 120051 ) and TrES-1 
ijCharbonneau et al.ll200E ) , that it has been possible to mea- 
sure the flux from these planets (for a review of t ransit- 
ing exoplanet science see ICharbonneau et alj [2007). The 
depth of the secondary eclipse, however, only characterizes 
the day-side flux from the planet. To fully characterize the 
planet's longitudinal temperature profile, observations must 
be made at a variety of plan etary phases. Th e first such 
observations w e re rep orted by lHarrington et all (|2006h and 
iKnutson et al.l (|2007l ). In this paper we present similar ob- 
servations of three other planetary systems. We briefly re- 
view models of Hot Jupiter phase variations in §1.2 and §1.3, 
and the observational constraints in §1.4. In §2 we present 
our observations and in §3 we explain how they constrain 
the properties of Hot Jupiters. We state our conclusions in 
§4. 
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1.2 Toy Model 

For the purposes of this paper, the phase function of a plane- 
tary system is taken to be the relative flux contribution from 
the planet as a function of its orbital position. Here we de- 
scribe a toy model of Hot Jupiter phase functions, based on 
a planet in a circular orbit with two isothermal hemispheres 
and neglecting the planet's remnant heat of formation, an 
effect ive temperature of Tf orm m 50-75 K (|Burrows et al.l 
2006). The Bond albedo, Ab, is the fraction of the incident 
stellar flux that is reflected by the planet. A planet's energy 
budget is largely determined by its Bond albedo. Hence, if 
we treat the host star as a blackbody, the equilibrium tem- 
perature of a planet's day and night faces is, respectively 
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where R* and T e s are the star's radius and effective tem- 
perature, a is the planet's semi-major axis and P n quantifies 
the portion of the absorbed stellar energy advected to the 
planet's night-side (P„ = for no redistrib ution; P n = 0.5 
for full redistribution) (|Burrows et al.ll2006l ). 

The depth of secondary eclipse (for edge-on systems) is 
approximately the ratio of the planet's day-s ide flux, -Fdav, 
to the stellar flux, F* . In light of the models of lFortnev et al.l 
(2006), we assume that the planet's day and night-side ra- 
diate as blackbodies [F — B\(T)) and that the host star 
emits with a mid-IR brightness temperature Tbright (we use 
Tbrigh t = 0.87T crl at 8 /im based on the models of iKuruczl 
1 19791 ). The resulting planet-star flux ratio is: 
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where R p is the planetary radius and the ratios R p /a and 
Rp/R* can be well-constrained if the planet transits its host 
star. The wavelength-dependent geometric albedo, A g , is the 
ratio of the total brightness of the planet at full phase to 
that of a Lambertian disk with the same cross-section. The 
second term (thermal emission from the planet) can be ne- 
glected at optical wavelengths. 

In the infrared, we may neglect the contribution from re- 
flected starlight. The peak-to-trough amplitude of the phase 
function for a transiting systems then becomes: 

Jday — Jnight _ BxjTday) ~ Bx{T n j ght) (Rv 
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where F n i g ht is the flux from the planet's night-side and 
where we have neglected the planet's flux contribution (< 
1%) in the denominator. Note that the absolute flux from 
the planet is a minimum during the secondary eclipse, ir- 
respective of phase variations, since the planet is entirely 
hidden from view at that time. The normalized phase vari- 
ations, AF(<f>) I (F) , for a system with inclination i is: 
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where the <j> is the planet's orbital phase (<f> — at inferior 
conjunction; <j> = ir at superior conjunction) and (F) is the 
time-averaged flux of the system. 



Measurements of secondary eclipse depth thus constrain 
a combination of Ab, A g and P n in this toy model. Ob- 
servations of the phase function amplitude constrain these 
same parameters, as well as sin i and R p in the case of non- 
transiting systems. 



1.3 Numerical Models 

In addition to the toy model outlined above, there have 
been numerous more sophisticated attempts to model the 
atmospheres of Hot Jupiters and thus to predict the depth 
of secondary eclipses as well as the full phase function of 
specific extrasolar systems. We summarize the most per- 
tinent aspects o f these models here ( for a more thorough 
background see iBurrows et al. I l2006h . ICooper fc Showmanl 
(2005], l2006h determine energy recirculation and predict 
global temperature contrasts of ~500 K at the photosphere. 
They also predict a noticeable offset betw een the phase 
functi on peak and the time of opposition. iFortnev et al.l 
|2005l ) construct models with either no or full recircula- 
tion and find that t he latter is a better fit to secondary 
eclipse obs ervations. [ Burrows et al.l l|2005l ). iBarman et al.l 
|2005l ) and ISeager et all ([20051 7 also run simulations with 
different amounts of recirculation and find that signifi- 
cant redistributio n of heat mu s t be occurring to explain 
the observations. IFortnev et al] ([2006T I run models to cal- 
culate emitted spectra as a function of orbital phase for 
both equilibrium and non-equilibrium chemistry. They find 
that non-equilibrium chemistry leads to significantly smaller 
phase variations and that current observations can be ex- 
plained with a black body-like planetary spectrum. Finally, 
IBurrows et al.l l|2006h use the factor P n to characterize the 
efficiency of energy recirculation and construct models for 
five different values of the parameter. They find that cur- 
rent observations rule out models with no recirculation and 
are in fact consistent with full recirculation (although they 
admit that current observations are not precise enough to 
determine the exact degree of recirculation). 



1.4 Observational Constraints 

Secondary eclipses have been reported for four short-period 
extrasolar planetary systems and a phase function has been 
reported in a different, n on-transiting syst e m. W e sum- 
marize these results below. iRichardson et al] (|2003l ) placed 
an upper limit of 3 x 10~ 4 on the relative secondary 
eclipse depth at 2.2 fim using the SpeX instrument at the 
NASA Infrared Telescope Facility. The secondary eclipse of 
HD 209458, observed with the Multiband Imaging Photome- 
ter for Spitzer (MI PS) yielded a relativ e depth at 24 fim of 
0.0026 ± 0.00044 (|Deming et al.ll2005l ). Similarly, the sec- 
ondary eclipse of TrES-1 was observed in the IRAC 4.5 and 
8 fim bands, yielding eclipse dept hs of .00066 ±0 .00013 an d 
0.00225 ± 0.00036, respectively j|Charbonneau et al l [20051 ). 
The latest syst em to have been dete cted in secondary eclipse 
is HD 189733 ([Deming et al.ll2006h . Observations with the 
Spitzer Infrared Spectrograph (IRS) peak-up imager at 
16 fim give an eclipse depth of 0.00551 ± 0.00030. In ad- 
dition to infrared detections, the non-detection of a sec- 
ondary eclipse for HD 209458 at optical wavebands has con- 
strained the 500 nm secondary eclipse depth to be less than 
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Figure 1. Sample 8 /im observation of science target, 51 Peg, 
and flux standard, HD 217636, showing the flux ramp-up. This 
observation began at MJD = 53728.447, or <j> = 1.8.87T. Each point 
represents 64 exposures, each of 0.1 sec. The error bars, which 
have an average value of 3.3 X 10 3 , correspond to the standard 
deviation in the counts from the 64 images. For comparison, one 
would expect uncertainties of roughly 2.2 X 10 3 in the Poisson 
noise regime. 

4.88 x 10~ 5 at the Iff le vel and less than 1.34 x 10~ 4 at the So- 
level (jRowe et al 1l2006l ). Ground-based K-band (2.2 /im) ob- 
servations of OGLE-TR-113b tent atively show a secondar y 
eclipse depth of 0.0017 ± 0.0005 ([Snellen fc Covinoll2007l ). 
MIPS 24 /im imaging of v Andromeda has yielded a phase 
function with peak-to-trough amplitude o f 0.0029 ± 0.0007 
and w eak evidence of a small phase offset (jHarrington et all 
2006). This surprisingly large amplitude suggests that the 
planet has parameters Ab ~ 0, P„ » and i ^ 30°. Con- 
tinuous IRAC 8 /xm monitoring of HD 189733 for half of a 
planetary orbit led to the detection of a phas e function with 
a pe ak-to-trough amplitude of 0.00121(4) jKnutson et al.l 
120071 ). This modest phase function, roughly one third the rel- 
ative depth of secondary eclipse, indicates that much of the 
incident ste llar energy is be i ng rec irculated to the planet's 
night side. iGrillmair et al.l |2006) used IRS to obtain a 
mid-IR spectrum for HD 189733b, which showed no evi- 
dence for absorpt i on ba nds due to either water or methane. 
I Richardson et al.l (|2007T J obtained the relative mid-IR spec- 
trum of HD 209458b and similarly concluded that such ab- 
sorption features could be excluded; they also identified can- 
didate emission features at two wavelengths. 



2 DESCRIPTION OF OBSERVATIONS 
2.1 Photometry 

We acquired Spitzer photometrjQ with the IRAC 3.6 (im 
or 4.5 /xm, and 8 /im detectors in sub-array mode 
(32x32 pixels, or 38.3x38.3 arcseconds) to avoid saturation 
( Fazi o et al ] |2004l ). Observations were made at eight differ- 
ent phases for each planetary system with follow-on time 
constraints. Each observation consists of four nods to the 
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target star interleaved with four nods to a nearby flux stan- 
dard of comparable infrared brightness. The number and du- 
ration of exposures for each of our targets are summarized 
in Table [1] We convert fluxes from the Spitzer units of spe- 
cific intensity (MJy/sr) to photon counts, and use aperture 
photometry to recover the counts for the science and stan- 
dard stars. We determine the PSF centroid using a centre- 
of-light function, then use the IDL routine APER to extract 
the source flux. After testing many combinations of aperture 
and sky annulus radii for robustness, we settle on a 7 pixel 
aperture radius and a sky annulus centred on the target star, 
with an inner radius of 10 pixels and an outer radius of 15 
pixels. 



2.2 Calibration 

A single observation lasts only 20 minutes, compared to 
planetary orbital periods of days, so the flux is not ex- 
pected to vary noticeably within a single observation. Un- 
fortunately, the measured signal in the 3.6 /im and 4.5 /im 
bands exhibits large (few %) fluctuations. The pointing 
drifts between nods are sub-pixel yet other researchers have 
found that the IRAC 4.5 /im In Sb array shows significant 
intra-pixel sensitivity variati ons (jCharbonneau et al.ll2005l ; 
iMorales-Calderon et al] [2006). The 3.6 /xm detector, also an 
InSb array, likely suffers from the same problem. We do 
not have sufficient images to properly map out the position- 
sensitivity relation, and this severely limits the quality of 
our 3.6 /im and 4.5 /im data. The 8 /im observations of 
51 Peg and HD 179949 show a smooth increase in gain, 
with a similar shape for both the science target and flux 
standard, indicating an instrumental effect, as shown in Fig- 
ure [1] HD 209458, a somewhat fainter source, does not ex- 
hibit such a marked ramp-up. The illumination-dependent 
n ature of this instrumen tal effect, described in greater detail 
in lKnutson et all (|2007l ). is consistent with charge trapping, 
in that pixel sensitivities are initially deficient but asymptote 
to a constant value as the pixels are illuminated repeatedly. 

For each observation, we determine the total counts by 
multiplying the median count by the number of frames, thus 
ignoring the minority of frames affected by cosmic rays. We 
repeat this median averaging for all the exposures within 
each nod (with duration T no d) and take the standard de- 
viation in the counts from each exposure multiplied by the 
square-root of the number of exposures per nod to be the un- 
certainty in the total counts. We find these uncertainties to 
be 4-5 times larger than the expectation from Poisson noise 
at 3.6 /im and 4.5 /im, and roughly twice Poisson noise at 
8 /im. At the shorter wavebands, as well as for the 8 /im 
observations of HD 209458, the points from different nods 
show considerable scatter but no obvious trend, so we add 
up the counts from all the nods separately for the science 
and reference star, then take the ratio of the two sums to 
determine the relative flux at that epoch. Noticing that the 
ramp-up at 8 /im for 51 Peg and HD 179949 levels off to- 
wards the end of any given observation, we take the ratio 
of the summed counts for the two final nods as the relative 
flux at that epoch. In all cases we compute the uncertainty 
in the flux ratio by taking the quadratic sum of the rela- 
tive uncertainties in the science and standard counts. The 
average flux measured for each star can be found in Table Q] 
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Table 1. Number and Duration of Exposures and Total Time on Target, Mid-IR Fluxes and 
Phase Variability of Target Stars 



Star 


A [|im] 


T CX p [sec] 


iVexp 


r nod [sec] 


F*(cr) [mjy] 


AF/ <F> 


51 Peg 


4.5 


0.02 


1728 


34.56 


533.1(9) 


< 0.017 




8 


0.1 


832 


83.2 


186.52(8) 


< 0.0007 


HD 217636 


4.5 


0.02 


1728 


34.56 


517.6(9) 


N/A 




8 


0.1 


832 


83.2 


198.2(1) 


N/A 


HD 179949 


3.6 


0.02 


2560 


51.2 


331(11) 


< 0.019 




8 


0.4 


576 


230.4 


73.72(4) 


0.00141(33) 


SAO 187890 


3.6 


0.02 


2560 


51.2 


367(2) 


N/A 




8 


0.1 


576 


230.4 


83.51(5) 


N/A 


HD 209458 


3.6 


0.1 


320 


32 


92.2(3) 


< 0.015 




8 


0.1 


320 


32 


20.75(4) 


< 0.0015 


HD 210483 


3.6 


0.1 


320 


32 


114.1(2) 


N/A 




8 


0.1 


320 


32 


25.49(5) 


N/A 



The wavelength of the observations is A, r exp is the duration of a single exposure, N exp is the 
number of exposures per nod, T no£ ; is the total time on target per nod, F*(<r) is the average 
flux of the system and the associated error, while AF/F r is the 2a upper limit, or the best-fit 
peak-to-trough amplitude and associated la uncertainty for phase variations, as appropriate. 
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Figure 2. The two bottom panels show the 8 /im light curves 
for the science target, 51 Peg, and the flux standard, HD 217636. 
The top panel shows the relative light curve (the ratio of the 
two lower plots). The phase of the planet (0 from §1.2) is at 
inferior conjunction, and 7r at superior conjunction. The la error 
bars have an average size of 0.00029 and represent the standard 
deviation in the mean of the binned counts. The curve is the best- 
fit sinusoid (with fixed period and no phase offset), as determined 
by a 10 5 -step Markov Chain Monte Carlo. 



Figure 3. The two bottom panels show the 8 fim light curves for 
the science target, HD 179949, and the flux standard SAO 187890. 
The top panel shows the relative light curve at for HD 179949. 
The phase of the planet (<j> from §1.2) is at inferior conjunction, 
and 7r at superior conjunction. The la error bars have an average 
size of 0.00028 and represent the standard deviation in the mean 
of the binned counts. The curve is the best-fit sinusoid (with fixed 
period and no phase offset), as determined by a 10 5 -step Markov 
Chain Monte Carlo. 



3 LIGHT CURVES & IMPLICATIONS 

We use ephemerides from lButler et ail {2006 ), for which pro- 
jected errors at the epochs of our observations are less than 
60 min. for all three systems. We neglect the correction for 
the light travel time across the Solar System, since the size 
of this correction (less than 8 minutes) is less than the uncer- 
tainty in the orbital phase. By matching the planet phases 
with the relative flux at that epoch, we construct a rela- 
tive light curve for each system, shown in the top panel of 
Figures [2}[4] (the 3.6 /im and 4.5 light curves are not 



plotted due to prohibitively large systematic uncertainties). 
HD 179949 and 51 Peg are non-transiting systems, while the 
observations of HD 209458 were taken outside of transit or 
secondary eclipse (represented by shaded areas in Figure [4}. 

We use a 10 -step Marko v Chain M onte Carlo with 
Metropolis-Hastings algorithm (tForc feOoBT) to determine the 
best-fit sinusoid, which is shown as the curve in the top panel 
of Figures [2}[4] The amplitude and constant offset of the si- 
nusoid are allowed to vary, while the period is fixed at 27r 
and the phase offset is set to zero, in ac c ordan ce with §1.2, 
as well as the results of lHarrington et al.l (|2006l ) . None of the 
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Figure 4. The two bottom panels show the 8 /jm light curves 
for the science target, HD 209458, and the flux standard star, 
HD 210483. The top panel shows the relative light curve. The 
phase of the planet (0 from §1.2) is at inferior conjunction, and 
7r at superior conjunction. The greyed-out areas correspond to a 
transit or secondary eclipse. The la error bars have an average 
size of 0.00070 and represent the standard deviation in the mean 
of the binned counts. The curve is the best-fit sinusoid (with fixed 
period and no phase offset), as determined by a 10 5 -step Markov 
Chain Monte Carlo. 

3.6 /im and 4.5 fim light curves, or indeed the 8 /im light 
curve for 51 Peg or HD 209458 show a variation that is better 
fit by a sinusoid than a straight line. We do, however, detect 
an 8 fim phase function with peak-to-trough amplitude of 
0.00141(33) for HD 179949, compared to an amplitude of 
0.0031 we would expect for an edge-on planet with a radius 
of 1.2Rj with an albedo of and a recirculation parameter of 
P n =0. For the two systems where we do not detect a phase 
function, the Markov Chain produces the cumulative distri- 
bution functions shown in Figures [5] and [6] These plots show 
the full range of theoretically allowed peak-to-trough phase 
function amplitudes on the x-axis, while the y-axis shows 
the probability that the system has an amplitude less than 
a given amplitude. The 1, 2 and 3<r upper limits are marked 
for convenience. The 2a upper-limits on phase variations at 
8 fim are 0.0007 for 51 Peg (compared to an amplitude of 
0.0026 for an edge-on planet with R — 1.2Rj with an albedo 
of and a recirculation parameter of P n = 0), and 0.0015 
for HD 209458 (compared to an amplitude of 0.0030 if the 
planet has an albedo of and a recirculation parameter of 
P n =0). Our complete results are summarized in Table [1] 

We recognize that v Andromeda and HD 179949 are the 
only t wo planet-bearing systems studied by IShkolnik et all 
(2005) to exhibit chromospheric activity in phase with the 
planetary orbit. It is therefore possible that the mid-IR 
phase variations that have been observed in these systems 
are a consequence of stellar activity and not the diurnal flux 
contrast of the planet. We do not believe that chromospheric 
activity is responsible for the phase function we observe in 
HP 179949, however, as the cycle in chromospheric activity 
IShkolnik et alj l|2005l ) observed in this system was 60° out 
of phase with the planetary orbit, and this pattern persisted 
for at least the year they observed it. 

If the phase variation of HD 179949 is taken at face 
value, it constrains the properties of this system, as ex- 
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Figure 5. Cumulative distribution function for the peak-to- 
trough amplitude of the phase function of 51 Peg based on Markov 
Chain Monte Carlo analysis of IRAC 8 fim images. The x-axis 
spans the full range of theoretically allowed amplitudes for this 
system. 

plained in Section 11.21 Constraining the planetary radius 
to within Rj < R p < 1.2Rj, the observed phase variation 
places a la upper-limit of P„ < 0.21. Further constraining 
the albedo t o Ab = 0, c onsist ent with the value of A g — 
0.03(5) that iRowe et al.l (|2007r i measured for HD 209458b, 
we obtain the exclusion plot shown in Figure [7] Our model 
thus places a la limit of M < 2AMj on the mass of 
HD 179949b. 

Because of the edge-on configuration of HD 209458, 
our upper-limit can be com bined with the c onstr aint on 
the planet/star radius ratio of iMandel fc Ago] (|2002l ). yield- 
ing constraints on the albedo and recirculation of the Hot 
Jupiter's atmosphere (s.f. § ll.2|l . In particular, the Ab — 0, 
P n — configuration of our toy model is excluded at the 3a 
level, as shown in Figure(8] More importantly, we exclude the 
much more realistic atmospheric models of iBurrows et al.l 
(2006) for P n ^ 0.4 by more than 2a, since their models 
predict stronger phase variations than our toy model for the 
same P n . 



4 CONCLUSIONS 

We have acquired Spitzer IRAC 3.6 /im or 4.5 /im, and 
8 fim images of three known Hot Jupiter systems. The use of 
standard stars has allowed us to achieve better than 0.05% 
relative photometry with multi-epoch IRAC 8 fim observa- 
tions. Intra-pixel sensitivity variations limit our photome- 
try at 3.6 /xm and 4.5 /im. We detect a 8 /im phase func- 
tion for HD 179949 and determine upper limits for the am- 
plitude of phase variations in the two other observed sys- 
tems. We analyze these light curves using a toy model of 
the planets and constrain the properties of the planets. 
Our observations suggest that Hot Jupiters do not all ex- 
hibit the large mid-infra red diurnal flux contrast detected by 
lHarrington et ail (|2006l '). Even within our modest sample of 
three planets, it appears that HD 179949b and HD 209458b 
exhibit rather different degrees of recirculation. The data 
for HD 179949 and v Andromeda provide the most un- 
certain constraints, however, while the two transiting sys- 
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Figure 6. Cumulative distribution function for the peak-to- 
trough amplitude of the phase function of HD 209458 based on 
Markov Chain Monte Carlo analysis of IRAC 8 /an images. The 
x-axis spans the full range of theoretically allowed amplitudes for 
this system. 
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Figure 7. Exclusion plot for the mass, M/Mj, and recirculation, 
P n , of HD 179949b based on its 8 /im phase function. P n = 
corresponds to no redistribution, while P n = 0.5 corresponds 
to full redistribution. Lower masses correspond to inclinations 
near it/ 2, while greater masses correspond to more face-on con- 
figurations. The radius of planet is assumed to be in the range 
Rj < R p < 1.2Rj, and its Bond Albedo is set to 0. The la and 
2a confidence intervals for large phase function amplitudes lie on 
top of each other so the former is omitted for clarity. 

tems with w ell-constrained phase functions (HD 209458 and 
HD 189733, iKnutson et~aDl2007r i both point toward a large 
value of P n > 0.35, so it remains to be seen if the obser- 
vations of the non-transiting planets hold up under further 
scrutiny. 
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Figure 8. Exclusion plot for the Bond albedo, Ab, and recircula- 
tion, P n of HD 209458b based on its 8 /jm phase function. P n = 
corresponds to no redistribution, while P n = 0.5 corresponds to 
full redistribution. Solid lines correspond to upper limits in this 
work. The dotted lines corresponds to constraints from the sec- 
ondary eclipse measurement (D. Charbonnea u, personal c o mmu - 
nication) , while the dashed line is from the iRowe et al.l d2007f ) 
albedo upper limit. The bottom-left corner of the plot (the no 
albedo, no recirculation configuration) is excluded at the 3c level 
by our data, while the shaded region in the top left of the figure 
is within all the la contours. 
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